1. Field of the Invention
The present invention relates to a method for updating a geological reservoir model representative of the structure and of the behavior of a petroleum reservoir, by integration of dynamic data.
More particularly, the invention applies to the construction of a map associated with the geological reservoir model, representative of the static and dynamic petrophysical properties and of their spatial variabilities provides engineers a means allowing better estimation of the reserves of a petroleum reservoir and to optimize recovery by selecting a suitable production scheme.
2. Description of the Prior Art
Optimization and development of petroleum reservoirs is based on the most accurate possible description of the structure and of the behavior of the reservoir A tool is used which accounts for these two aspects in an approximate way: a geological reservoir model. A geological reservoir model best accounts for the structure and the behavior of the reservoir. The model includes a grid pattern that forms the frame of the reservoir and that has to be representative of the structure, and petrophysical property maps associated with this grid, that have to be representative of the behavior. This association assigns a petrophysical value obtained from maps to each grid cell of the model.
A petrophysical property map describes the spatial distribution or spatial structure, in an underground zone, of continuous or discrete petrophysical properties such as permeability, porosity, lithology type, . . . Conventionally, these maps are obtained from stochastic methods and are then referred to as stochastic models. In a stochastic context, the expression realization is used rather than numerical model.
The quality of the optimization and of the development of a petroleum reservoir thus directly depends on the accuracy of the stochastic realizations (petrophysical property maps). It is therefore advisable to work out stochastic realizations and therefore, more generally, reservoir models that are as coherent as possible with all the data collected (well, seismic, laboratory data, . . . ).
The data available for constraining stochastic realizations are referred to as static or dynamic. A datum is static if it corresponds to a measurement of the property modelled at a given point and if it is independent of time. Permeability measurements performed in the laboratory on rock samples or logs measured along wells are static data. A datum is dynamic if it depends on time, it is linked with the property modelled without being a direct measurement thereof. Production data and 4D seismic data, that vary with fluid flows, are dynamic data. Since the data are insufficient to allow deterministic description of the spatial distribution of the property considered, stochastic modelling techniques are most often based on geostatistical techniques, which provide a family of numerical stochastic models associated with the geological reservoir model and referred to as realizations.
In a stochastic context, data describing the geology of the medium define a random function. For a single random function, there is an infinity of possible realizations. All these realizations are however not compatible with dynamic data.
Static and dynamic data are not integrated in the same way in the stochastic realization. Integration of static data is carried out upon generation of the realization whereas integration of dynamic data goes through the solution of an inverse problem involving a flow simulator. Dealing with the inverse problem first requires defining a function referred to as objective function, or cost function, that measures the relevance of the realization or of the reservoir model proposed. In the first research work devoted to this subject, the objective function was a direct measurement of the difference between the dynamic data collected in the field and the corresponding responses obtained by simulation:
      J    ⁡          (      y      )        =            ∑      i        ⁢                            w          i                ⁡                  (                                    d                              sim                i                                      -                          d                              obs                i                                              )                    2      y is the realization considered and J(y) the value of the objective function for this realization. The w are weighting coefficients. The dobs are the dynamic data collected and the dsim are the simulated corresponding responses. The quantity g is the operator that goes from the space of the non-constrained realizations to the space of the dynamic data: dsim=g(y). Minimizing this objective function leads to determination of a realization y that reproduces all of the dynamic data as well as possible. Unfortunately, the spatial structure of realization y thus obtained is generally no longer coherent with that of the initial realization, that is with the geologic data.
A framework more suited to the definition of the objective function is provided by the Bayesian approach. A priori information is then added to the objective function. This approach is described in the document as follows:                Tarantola, A., 1987, “ Inverse Problem Theory—Methods for Data Fitting and Model Parameter Estimation”., Elsevier, Amsterdam, 613 p.        
The objective function is then expressed as follows:
      J    ⁡          (      y      )        =                    1        2            ⁢                        (                                    g              ⁡                              (                y                )                                      -                          d              obs                                )                t            ⁢                        C          D                      -            1                          ⁡                  (                                    g              ⁡                              (                y                )                                      -                          d              obs                                )                      +                  1        2            ⁢                        (                      y            -                          y              o                                )                t            ⁢                        C          Y                      -            1                          ⁡                  (                      y            -                          y              o                                )                    
The first term of the objective function deals with the likelihood constraint: it measures the difference between the dynamic data observed in the field and the equivalent data obtained by numerical simulation. The second term corresponds to the a priori constraint: it quantifies the difference between the a priori reservoir model yo, deduced from the a priori geologic information, and the proposed reservoir model y. The covariance matrix CD characterizes the experimental and theoretical uncertainties, whereas Cy concerns the uncertainty on the a priori model. Minimization of the objective function then provides a model y that is as close as possible to the a priori model and such that data dsim simulated for this model are close to the data measured in the field.
Calculation of the objective function can be difficult, notably on account of the two covariance matrices that have to be inverted. In general, the covariance matrix relative to the data is assumed to be diagonal and it is readily inverted. This hypothesis is most certainly questionable depending on the case considered, for example with the 4D seismic method, but it is not challenged here. The case of the covariance matrix relative to model y is more awkward. In fact, the dimension of a priori covariance matrix Cy is the length of the vector y and calculation of its inverse is often impossible for models comprising a very large number of grid cells. To date, to simplify taking account of the a priori constraint, essentially three types of approach have been developed. The first type of approach is based on the division of covariance matrix Cy into subspaces: the least influent components are disregarded, which allows the number of variables to be reduced. This technique is described in the document as follows:    Reynolds, A. C., He, N., Chu, L., and Oliver, D. S., 1995, “Reparameterization Techniques for Generating Reservoir Descriptions Conditioned to Variograms and Well-Test Pressure Data.”, SPE Annual Technical Conference and Exhibition, Dallas, Tex., 22-25 October, SPE 30588, p. 609 -624.
The second type of approach is based on the mathematical modelling of errors in the parameters space by a Gaussian random function of mean zero and of exponential covariance along preferred directions. There are thus mathematical properties allowing to analytically calculate the inverse of the a priori covariance matrix. This method is described, for example, in French Patent 2,635,197 and corresponding U.S. Pat. No. 4,972,383.
Finally, the third type of approach involves geostatistical parameterization, with in particular the pilot point method introduced in the document as follows:    Marsily, G. de, 1978, “De I'identification desSystémes Hydrologiques”., Thése de Doctorat d'Etat, Université Paris VI, France.and the gradual deformation method described in the document as follows:    Hu, L. Y., 2000, “Gradual Deformation and Iterative Calibration of Gaussian-Related Stochastic Models.”, Math. Geology, v. 32, no. 1, p. 87-108.
The pilot point method was first introduced within the context of estimation prior to being extended to the conditioning of stochastic realizations by dynamic data by:    RamaRao, B. S., LaVenue, A. M., de Marsilly, G., and Marietta, M. G., 1995, “Pilot Point Methodology for Automated Calibration of an Ensemble of Conditionally Simulated Transmissivity Fields. 1. Theory and Computational Experiments.”, Water Res. Res., 31(3), 475-493.
This technique allows local deformation of realizations from a reduced number of parameters while respecting the spatial variability of the property modelled (permeability, porosity, velocity, . . . ). In short, in order to modify the realization, a set of points (or cells) referred to as pilot points, whose values can be changed, is selected. The perturbation generated at these points spreads over the entire realization by kriging according to the expression:yx(x)=ydK(x)+[y(x)−yk(x)]y is a non-constrained realization. The quantity ydK results from kriging of the static data and the values of the pilot points, and yK from kriging of the values of y at the same points. The quantity yc is a constrained realization that honors the spatial variability model and the values of the static data, as well as the values selected for the pilot points. In other words, the pilot points are regarded as static data that are used to constrain the realization. The values of the pilot points are the parameters of the optimization problem. These “pseudo”-data, unlike the real static data, are not stationary: they can vary during the optimization procedure so as to reduce the objective function. Since the modifications are propagated to the entire realization by kriging, conservation of the spatial variability model is ensured. Adding the term relative to the a priori constraint in the objective function is therefore considered to be redundant. This objective function then comes down to the single term measuring the difference between the real dynamic data and the corresponding responses obtained by simulation. It is no longer necessary to determine the inverse of the a priori covariance matrix. This property is fundamental and characteristic of the pilot point method.
However, the pilot point method can involve numerical artifacts. In some cases, minimization of the objective function goes through the assignment of extreme values to the pilot points. We can then have either too high values or too low values, which no longer make sense, physically speaking. To avoid these artifacts, constraints by inequalities can be integrated in the optimization procedure: the variations of the parameters are then bounded. This technique is described more precisely in the aforementioned document by RamaRao et al. (1995). Besides, the values of the pilot points are adjusted independently of one another. The pilot points therefore have to be separated by a distance that is greater than or equal to the correlation length. If this minimum distance is not respected, the pilot point method does not guarantee preservation of the spatial variability model.
The gradual deformation method was initially proposed to continuously modify Gaussian stochastic processes (random processes). It is a geostatistical parameterization technique allowing deformation of a realization of a reservoir, model comprising any number of grid cells from a reduced number of parameters, while respecting the spatial variability model. The basic idea is that the sum of Gaussian random functions is a Gaussian random function.
The simplest gradual deformation scheme adds up two multi-Gaussian random functions. Let Y1 and Y2 be two such functions, independent and stationary of order 2. They are also assumed to have the same mean (yo), variance and covariance model. A new random function Y(t) is constructed by combining Y1 and Y2 according to the expression as follows:Y(t)−yo=└Y1−yo┘cos(t)+└Y2−yO┘sin(t)
It can be shown that, for any deformation coefficient t, Y has the same mean, variance and spatial variability model as Y1 and Y2. In fact, the sum of the coefficients squared, i.e. cos2(t)+sin2(t), is 1. According to this combination principle, construction, from two realizations independent of Y1 and Y2, denoted by y1 and y2, of a chain of realizations depending only on deformation parameter t is possible.y(t)−yo=└y1−yo┘cos(t)+└y2−yo┘sin(t)
This chain of realizations goes through y1 and y2. When t is 0, y is y1. When t is □/2, y is y2. By continuously varying the deformation coefficient t from 0, continuous deformation of realization y1 is simulated. By varying continuously deformation coefficient t from 0, continuous deformation of realization y1, taken as the initial realization, is simulated. An essential point is that, for any value of t, realization y is multi-Gaussian and respects the mean, variance and spatial variability model of y1 and y2.
When the gradual deformation method is integrated in an optimization process as a parameterization technique, the objective function to be minimized becomes:
      J    ⁡          (      t      )        =                              1          2                ⁡                  [                                    g              ⁡                              (                t                )                                      -                          d              obs                                ]                    t        ⁢                            C          D                      -            1                          ⁡                  [                                    g              ⁡                              (                t                )                                      -                          d              obs                                ]                    .      
In fact, as for the pilot point method, it seems redundant to add the a priori constraint in the objective function because parameterization intrinsically preserves the spatial variability model. Vector t combines the various deformation coefficients. These coefficients are the parameters of the optimization problem. These deformation parameters have to be identified so as to reduce the objective function as much as possible.
The gradual deformation method as presented here implies that the entire realization is deformed to minimize the objective function. However, the gradual deformation method can also be applied locally. In this case, instead of combining realizations with a given spatial variability model, the Gaussian white noises used to generate these structured realizations are combined. More precisely, when the deformation is to be located in a given zone, the components (random numbers) of the Gaussian white noise assigned to the cells included in the zone to be deformed are gradually deformed. This technique is described in the document as follows:    Le Ravalec, M., Noetinger, B., and Hu, L.-Y., 2000, “The FFT Moving Average (FFT-MA)Generator: An Efficient Numerical Method for Generating and Conditioning Gaussian Simulations.”, Math. Geol., 32(6), 701-723.If necessary, the gradual deformation method can be applied to an isolated component of the Gaussian white noise. It then tends towards the pilot point method. An important difference has to be noted: the gradual deformation method prevents the modified point from taking extreme values. Furthermore, if several points are modified according to the gradual deformation method, the spatial correlations between the deformed points are taken into account. However, the gradual deformation method is negligible when it is applied to points. An optimization procedure involving this type of parameterization will therefore probably be much less efficient than the pilot point method.